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Abstract 

Cold nuclear matter effects on Jlxj/ production in proton-nucleus and nucleus-nucleus collisions are evaluated taking into account the 
specific y/i^-production kinematics at the partonic level, the shadowing of the initial parton distributions and the absorption in the nuclear 
matter. We consider two different parton processes for the cc-pair production: one with collinear gluons and a recoiling gluon in the final 
state and the other with initial gluons carrying intrinsic transverse momentum. Our results are compared to RHIC observables. The smaller 
values of the nuclear modification factor /J^^i in the forward rapidity region (with respect to the mid rapidity region) are partially explained, 
therefore potentially reducing the need for recombination effects. 
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1. Introduction 

The charmonium production in hadron collisions is an 
on-going major subject of investigations, on both experi- 
mental and theoretical sides. It has been widely studied 
in pp collisions; our understanding was recently reviewed 
in [1, 2]. It may also be used as a tool to probe the 
medium produced in nucleus-nucleus (AB) collisions (for 
a recent review, see [3] along with some perspectives for 
the LHC [4]). This medium is expected to be in a de- 
confined state of QCD matter - such as the Quark Gluon 
Plasma (QGP) - at high enough temperatures and densi- 
ties. The J /if) production should be sensitive to the QGP 
formation, due to competing effects such as a color De- 
bye screening suppression [5] or the so-called recombina- 
tion mechanism [6]. Recent results on J/i// production are 
available from the PHENIX experiment at the BNL Rel- 
ativistic Heavy Ion Collider (RHIC). They show a signi- 
ficant suppression of the J /iff yield in AuAu collisions at 
yfsNN - 200 GeV [7] compared to the expected yield from 
pp measurements [8]. However, the interpretation relies on 
a good understanding and a proper subtraction of the Cold 
Nuclear Matter (CNM) effects, known to impact the J/tfr 
production in proton(deuteron)-nucleus (pA, dA) collisions 
where the deconfinement can not be reached. Indeed, the 
pA data [9] obtained at the SPS energies can be described 
by assuming the break-up of the pre-resonant cc pair due to 
multiple scattering along its way to escape the nuclear envi- 
ronment - the so-called nuclear absorption. PHENIX data 
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on dAu collisions [10] have also revealed that CNM effects 
play an essential role at RHIC energy. These effects a priori 
include shadowing, i.e. the modification of the paiton dis- 
tribution of a nucleon in a nucleus, and final-state nuclear 
absorption. 

As we shall show thereafter, the impact of gluon sha- 
dowing depends on the partonic process producing the 
cc and then the J/i//. So far, the studies of pro- 
duction [11, 12, 13] with gluon shadowing relied on the as- 
sumption that the cc pair was produced by the fusion of two 
gluons carrying some intrinsic transverse momentum kj-- 
The partonic process being a 2 — > 1 scattering, the sum 
of the gluon intrinsic transverse momentum is transferred 
to the cc pair, thus to the J/i// since the soft hadronisation 
process does not modify the kinematics. This corresponds 
to the picture of the Colour Evaporation Model (CEM) at 
LO (see [1] and references therein). 

In such approaches, the transverse momentum of the J /iff 
entirely comes from the intrinsic transverse momentum of 
the initial gluons. This seems acceptable for the I0W-P7 re- 
gion and the origin of this intrinsic Pj can be paralleled to 
the increase of (P^) when going from pp to pA and for in- 
creasing atomic number A. This is known as the Cronin 
effect: the increase of (Pj) is believed to come from a 
broadening of the intrinsic transverse momentum distribu- 
tion, resulting from the multiple scatterings experienced by 
the initial gluon from the proton as it goes through the target 
nucleus before the heavy-quark production [14]. 

However, such an effect is not sufficient to describe the 
Pt spectrum of quarkonia produced in hadron collisions [1]. 
Most of the transverse momentum should have an extrinsic 
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origin, i.e. the Pj would be balanced by the emission 
of a recoiling particle in the final state. The J /iff would then 
be produced by gluon fusion and with emission of a hard 
final-state gluon. For the production of ^5 1 states - like the 
J /if/, such a 2 — > 2 partonic process is anyhow mandatory 
to satisfy C-parity conservation. 

It is among our purposes here to investigate the influence 
of such an emission on the kinematics of the J/i// production 
in p{d)A and AA collisions. Indeed, for a given J/i// mo- 
mentum (thus for fixed y and Pj), the processes discussed 
above, i.e. g+g cc ^ J/tfr(+X) and g+g — > Jj^+g, will 
proceed on the average from gluons with diff'erent Bjorken- 
X. Therefore, they will be affected by different shadowing 
corrections. From now on, we will refer to the former sce- 
nario as the intrinsic scheme, and to the latter as the ex- 
trinsic scheme. In the following, we shall consider them as 
distinct approaches^. 

In practice, we shall study the kinematic regime at RHIC 
at yfsNN - 200 GeV. For the extrinsic scheme, we shall 
consider the partonic diff'erential cross section forg + g^ 
J/ilf + g given in [15] which satisfactorily describes the data 
obtained in pp collisions at RHIC down to Pj ~ (see 
Fig. 1). For the intrinsic scheme, we shall follow the stu- 
dies [11, 12, 13] based on (2 — > l)-like processes where 
the momentum of the particles denoted hy X in g + g 
cc — > J/iff{+X) is neglected . Using a probabilistic Glauber 
Monte-Carlo code, we interface these production processes 
with CNM effects, such as shadowing and nuclear absorp- 
tion, in order to get the J/i// production cross sections for 
pA and AA collisions. We shall finally compare our results 
with the experimental measurements presently available at 
RHIC. 



2. The Monte-Carlo framework for J/iJ/ production 

To describe the J/t// production in nucleus collisions, 
our Monte-Carlo framework is based on the probabilistic 
Glauber model, the nuclear density profiles being defined 
with the Woods-Saxon parameterisation for any nucleus 
A > 2 and the Hulthen wavefunction for the deuteron [17]. 
The nucleon-nucleon inelastic cross section at -<i/snn = 
200 GeV is taken to o-nn = 42 mb and the maximal nucleon 
density to po = 0.17 nucleons/fm"^ . For each event (for each 
AB collision) at a random impact parameter b, the Glauber 
Monte-Carlo model allows us to determine the number of 
nucleons in the path of each incoming nucleon, therefore 
allowing us to easily derive the number A^coii of nucleon- 
nucleon collisions and the total number Np^n of nucleons 
participating into the collision. 



^In the extrinsic scheme and for this first study, we shall neglect the 
small kinematical effects of the yield from the decay of Xc produced by 
2 — » 1 processes. Indeed, from the 25 % of the Xc feeddown [16], it is 
reasonable to suppose that only half of it effectively proceeds via a gg —> 
Xc which is allowed at LO in ag only fovxd, not for^^l- 
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(a) Pt distributions from the i-channel cut contribu- 
tions compared to PHENIX pp data [8] in the forward 
and central rapidity regions. 
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(b) Rapidity spectrum from s-channel cut contributions compared 
to PHENIX pp data, the ad hoc double-Gaussian fit [8] and the 
predictions from the PYTHIA event generator 



Figure 1; Pj and y specti-a in pp at 



200 GeV. 



2.1. The CNM effects 

To get the J/i// yield in pA and AA collisions, a 
shadowing-correction factor has to be applied to the J/i// 
yield obtained from the simple superposition of the equi- 
valent number of pp collisions. This shadowing factor can 
be expressed in terms of the ratios Rf of the nuclear Parton 
Distribution Functions (nPDF) in a nucleon of a nucleus A 
to the PDF in the free nucleon. In the following, we will 
consider the evolution model EKS98 [18] which is recog- 
nised to be a reasonable compromise between DS [19] and 
EPS08 [20] as regards the strength of the gluon antishado- 
wing for instance. It provides Rf at a given initial value 
of jjf - the factorisation scale - and takes into account their 
evolution through the DGLAP equations. The nuclear ratios 
of the PDFs are expressed by: 



Rf(x,PF) = 



ff(x,PF) 
Af^"'^'""\x,pF) 



fi = q^q^g ■ 



(1) 



Within EKS98, these nuclear ratios are parameterized at 
some initial scale //^q = 2.25 GeV^ which is assumed large 
enough for perturbative DGLAP evolution to be applied. 
They are evolved at LO from yu^^ up to //^ (< lO'* GeV^) 
and are valid for x > 10"^. The numerical parameterisa- 
tion of Rf(x,pf) is given for all parton flavours. Here, we 



restrain our study to gluons since, at high energy, y/i/r is es- 
sentially produced through gluon fusion [1]. Usually, the 
spatial dependence of R'^{x,Hf) is not given. However, as 
we shall see in the next section, it can be included in our 
approach. 

The second CNM effect that we are going to take into 
account concerns the nuclear absorption. In the frame- 
work of the probabilistic Glauber model, this effect refers 
to the probability for the pre-resonant cc pair to survive to 
the propagation through the nuclear medium and is usually 
parametrised by introducing an effective absorption cross 

section (Tabs. 

2.2. The differential J /i^ -production cross section in 
AB collisions 

Our Glauber Monte-Carlo framework is aimed at nume- 
rically evaluating the differential y/i/r-production cross sec- 
tion in nucleus collisions, by exploring the whole physical 
phase space with a proper weighting of each point in this 
phase space. Within this Glauber-based calculation, the dif- 
ferential cross section for the production (via gluon fusion) 
of a quarkonium with momentum {y, Pj), in nucleus colli- 
sions at impact parameter b and for a nucleon-nucleon CM 
energy of yfsNN can be represented by a generic integral. 
It takes two different forms depending on the kinematics of 
the partonic process responsible for the y/f/r production. 

The intrinsic scheme corresponds to a 2 — > 1 partonic 
process, with initial gluons carrying a non-zero intrinsic 
transverse momentum. Following [13], we do not neglect 
the value of the 7/t/r's Pj in this simplified kinematics. In 
this scheme, the measurement of the Jj^ momentum com- 
pletely fixes the longitudinal momentum fraction carried by 
the initial partons: 



= exp (±y) = xl2(y,PT), 



JSnn 



(2) 



with the transverse mass mj = -^M^ + Pj, M being the 
y/i/' mass. 

Therefore, we can write 



J Intr. 



drAdzAdzBT^ix°,rA,ZA,I^F) (3) 



dy dPj db 

Thxl 4, ZB,^lF) 4') SA(fA,ZA)SB{rB, Zb) 



with 



'a ['b -b- r^] and za [zb] are the transverse and lon- 
gitudinal spatial locations of the initial parton (a gluon 
here) in the nucleus A [E]; it carries the longitudinal 
momentum fraction Xj {x^ and can be found with the 
probability Tg U^g] at a scale jU/r; 

The nuclear absorption is taken into account 
through 5a (rA, Za) = exp (-A crabs j^J t/z Pa('*a,z)) 
[S b^Tb, Zb)], which stands for the survival probability 
for a cc produced at the point {rA,ZA) [(''b,Zb)] to pass 
through the projectile and the target unscathed; 



• cr^'^"' (xp Xj) is the partonic cross section for the pro- 
cess g -I- g — > cc — > J/iff(+X) and is a function of 
Pj and y through jCj and x^. As we will show in the 
section 2.3, it can be extracted from experimental data 
along with the PDFs. 

In the extrinsic scheme, we deal with a 2 — > 2 partonic 
process with coUinear initial gluons and we have 



AB^J/if/X 



-- J'dxidx2 J" 



dy dPj db 

Tg{X2, rB, ZB,f^F)2sPT 
SA(r,ZA)SB(rB,ZB) , 



drAdzAdzBTg (xi , rA,ZA,l^F) 

do-gg^jlf+g 



di 



-5{s -t-u- M ) 



(4) 



with 



S — Sf\[js[X]X2, 



- M - x\ ^Sfit^mje^ , U - M - 
X2^SMMmje The four-momentum conservation 
- represented by the 6 function - explicitly results 
in a more complex expression of X2 as a function 
of (xi,y,f7-): 



X2 = 



Xinij ys^N^ 



(5) 



Equivalently, a similar expression can be written for xi 
as a function of (x2,y, Pt)', 

• da-gg^j/^+g/di can be computed in a priori different 
approaches that correctly describes the pp data. For 
now, we use the one obtained in [15], but others can 
be interfaced with our code. See section 2.3 for further 
details. 

Now, concerning 'Fgixi,rA,ZA,t^F), we assume that 
it can be factorised in the nuclear density distribution 
PAirA,ZA), the shadowing modification factor 'Rg{rA,x,pf) 
and the usual gluon PDFs g{x; fj.f): 

Tg{xi,rA,ZA;iJ-F) ^ pA(rA,ZA)'J<g(rA,xi,pf) g(xi;pf) . 

A priori, the modifications of the nPDFs should depend 
on the parton position (r, z) in the nucleus. Such infor- 
mation is not experimentally available. So the centrality 
dependence is not encoded in the EKS98 parametrisation. 
However, some approaches provide an Ansatz for such a 
dependence. Assuming that the inhomogeneous shadowing 
is proportional to the path length [21, 12], then 



Wa,x,pf) = 1 + [Rgix,fiF) - m> 



J dzpA(rA,z) 
J dzpA(0,z) 



(6) 



where R'^(x, fif) is the ratio nPDF/PDF given by the EKS98 
parametrisation for the gluon (see Eq. (1)) and A'p, is a nor- 
malisation factor, determined such that 



dzAPA(rA,ZAWg(rA,x,pF) = RJx,pf) 



The integral over z in Eq. (6) includes all the material tra- 
versed by the incident nucleon. This amounts to consider 
the incident parton as coherently interacting with all the tar- 
get partons along its path length. 

In the following, we shall set the scale yU/r in equal 
to the renormalisation scale of the partonic process and 
take the usual choice jUs = mj, mj providing a typical scale 
of the partonic process. 

2.3. The partonic cross sections 

Within the intrinsic scheme, the measurement of the dif- 
ferential cross section in pp on Fig. 1 (a) and (b) directly 
provides us with values for the product of the PDFs and the 
partonic cross section in Eq. (3). Indeed, 



^ Intr. 



(7) 



dydPT 

In order to evaluate the integral Eq. (3), we randomly pick 
y and Pj out of their respective distributions obtained from 
fits to the experimental pp spectra dcr^pp /dy dPr ■ For the 
present study, we use the fits to the y ^ and Pj spectra mea- 
sured by PHENIX [8] in pp collisions at y/s^ = 200 GeV 
as inputs of the Monte-Carlo. The azimuthal-angle tp of Pt 
in the (P^, P,,) plane is also random and follows a flat dis- 
tribution within [0, In]. As discussed previously, the know- 
ledge of {y,PT) unequivocally fixes the other variables (xi, 
X2 and fif - nij) needed to compute the shadowing correc- 
tion factors. 

On the other hand, in the extrinsic scheme, information 
from the data alone - the y and Pj spectra - is not suf- 
ficient to determine x\ and X2. Indeed, the presence of a 
final-state gluon authorises much more freedom to choose 
{x\,X2) for a given set (y, fj-). Even if kinematics deter- 
mine the physical phase space, models are anyhow manda- 
tory to compute the proper weighting of each kinemati- 
cally allowed {x\,X2). This weight is simply the differential 
cross section at the partonic level times the gluon PDFs, 
i.e. g(xi,iif)g(x2,iJ.F)do-gg^j/^+g/dydPTdxidx2. In the 
present implementation of our code, we are able to use the 
partonic differential cross section computed from any theo- 
retical approach. For now, we use the one from [15] which 
takes into account the .s-channel cut contributions [22] to 
the basic Color Singlet Model (CSM) [23] and satisfacto- 
rily describes the data down to very low Pj, where the bulk 
of the cross section lies. As shown on Fig. 1, this approach'' 
gives a fairly good description of both the Pj- and y depen- 
dence of the pp data at RHIC. 

To evaluate the integral in the extrinsic scheme, we also 
randomly pick y and Pj, but out of the distributions com- 
puted with the cross section computed as in [15]. For a 
given set (j, fr), the set (xi,X2) is randomly chosen in its 
kinematically allowed range and follows the distribution 
g{xi , fif) g(x2,iUF) da-gg^ji^+gldy dPTdx\dx2. 



The double gaussian parametrisation as quoted in [8]. 
''The determination of the two parameters of this approach [15] has 
been improved by fitting RHIC pp data (a = 3.2 and k = 6.3 GeV). 



2.4. Picturing intrinsic vs extrinsic scheme 

Fig. 2 shows the physical phase space in the {x2,y) plane 
for both schemes in JAu collisions. At a fixed value of y, 
they give quite different distributions of the Bjorken-jc of 
the initial gluons. The momentum of the final-state gluon 
considered in the extrinsic scheme results in a larger X2- 
range, whereas the intrinsic scheme much heavily favours a 
tighter band at low X2 ■ 




y -2 0.6 




Figure 2: Normalised X2 distribution versus rapidity in the 
intrinsic and extrinsic scenarios for similar Monte-Carlo 
Jjij/ statistics in c/Au collisions at yfsNN - 200 GeV. 



3. Results 

In the following, we present our results for the Jlxjj nu- 
clear modification factor: 



R 



dN- 



AB 



{N,,n)dN-'J/' 



(8) 



dN-'/^idNp'll') is the yield observed in AB (pp) cofli- 
sions and (A'coii) is the average number of nucleon-nucleon 
collisions occurring in one AB collision. In the absence of 
nuclear effects, R^b should equal unity. 

3.1. dAu collisions 

PHENIX measurements of Rd au [10] provides with a 
means to size-up the CNM effects at play at RHIC energy. 
We shall compare the CNM effects obtained in the intrinsic 
and extrinsic schemes to these data. 

Fig. 3a shows Rjau versus y. Let us first focus on the 
curves without nuclear absorption. The enhancement due to 
antishadowing is shifted to more negative y in the intrinsic 
scheme^ compared to the extrinsic scheme. This is easily 
explained. We saw in section 2.4 that the typical values of 
X are increased in the extrinsic scheme for given y's. A less 



'The dot dashed curve is obtained for a fixed Pf = and the dashed 
one with open squares for /if = mr using the distribution of the pp data 
in the forward rapidity region [13]. The results are in practice identical 
if we consider the distiibution in another v-region. In the following, all 
the curves covering the entire rapidity range for the intrinsic scheme are 
obtained as this dashed curve. 




Fi gure 3: J jip nuclear modification factor in t/Au collisions at ■sfsf^ — 200 GeV. 
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(c) RcuCii in the extrinsic scheme. 



Fi gure 4; A^pait dependence of the J I if/ nuclear modification factor in CuCu and AuAu collisions at ■^~si^ — 200 GeV. 



negative value of y is therefore required to obtain a value 
of X2 producing the maximum amount of antishadowing as 
seen on Fig. 3a. 

However, as usually done, it is necessary to add the nu- 
clear absorption. In the intrinsic scheme, we have used the 
value of (Tabs - 2.8 mb from Ref [10], where it was ob- 
tained by fitting the data with a shadowing-correction cal- 
culation equivalent to neglecting Pj in the intrinsic scheme. 
A higher value of crabs is required in the extrinsic scheme. 
We used (Tabs = 4.2 mb which gives a good agreement with 
PHENIX dAu data^. The resulting curves are also plotted 
on Fig. 3a. From now on, we will keep these values for the 
value of (Tabs in the respective schemes. 

Fig. 3b shows RdAu as a function of A^coii in the three rapi- 
dity windows. Both nuclear shadowing and absorption are 
included. Both schemes agree well with the data, especially 
at mid-y 



*A complete fitting procedure taking into account the experimental er- 
rors and their correlations is beyond the scope of this first analysis of the 
extrinsic effects. By limiting our analysis to a constant value for o"abs, we 
also disregard possible formation-time effects. 



Fig. 3c shows Rdku as a function of Pj in the three rapi- 
dity windows. The intrinsic and extrinsic approaches give 
similar results. They both give a reasonable agreement with 
the data. 

3.2. Nucleus-nucleus collisions 

We now extend our study to AA collisions, where the 
Jjij/ yield may be further affected by dense matter effects. 
In AuAu collisions [7], PHENIX has measured an unex- 
pected stronger y/i/r suppression at forward-y than at mid- 
y. A possible explanation lies in the recombination scena- 
rios [6]. In the following, we shall rather investigate if part 
of this rapidity-dependent suppression may be due to CNM 
effects. We shall also compare our results to PHENIX CuCu 
data [24]. They cover with better precision the lower A^pan 
range, where a smaller amount of dense matter effects is 
expected. 

Fig. 4a shows ^auAu versus A^pan in the intrinsic scheme, 
while Fig. 4b and Fig. 4c show i^cuCu and /?auAu for the 
extrinsic scheme. As regards /?cuCu and 7? auAu versus y, they 
are displayed on Fig. 5 for four different centrality bins. 



In the intrinsic case, /?auAu is nearly independent of y 
(Fig. 5 up). Indeed, for any value of y in the rapidity range 
-2 < y < 2, the suppression at y > in one nucleus is com- 
pensated by an enhancement at y < in the other nucleus, 
thus giving the same result as at mid-rapidity. 

On the other hand, in the extrinsic scenario, this cancel- 
lation is not as effective and both 7?cuCu and RauAu show 
a maximum at y = (Fig. 5). This is explicit in Fig. 4b 
and Fig. 4c where the curves for the central and forward 
rapidity ranges are shifted from each other, as the data are. 
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Figure 5: J/i// nuclear modification factor versus rapidity in 
AuAu (top) and CuCu (bottom) collisions. 



of the production at the partonic level. We have studied 
both intrinsic (2 1-like process) and extrinsic (2 — > 2- 
like process) production schemes. 

We have obtained different gluon-shadowing-induced ef- 
fects in dAu depending on the considered scheme. We have 
then included the nuclear absorption cross section needed 
to reproduce PHENIX c/Au data [10]. In the simplified 
kinematics of a 2 — > 1 process (previously considered 
in [11, 12, 13]), we have used cr^bs - 2.8 mb according 
to [10]. Within the extrinsic scheme, a larger break-up cross 
section is needed and we used cr^\,s - 4.2 mb. 

Concerning nucleus-nucleus collisions, we have studied 
both CuCu and AuAu collisions in order to compare our 
prediction with PHENIX measurements [7, 24]. Within 
the extrinsic scheme, we observed a rapidity dependence 
of ^cuCu and /?auAu, in the same direction as the one seen in 
the data. 

In the near future, we plan to extend our investigations to 
the LHC energies, where we could consider the production 
of T in pA and AA collisions at nonzero Pj by interfa- 
cing partonic matrix elements obtained at NLO [25, 26] 
and NNLO*[27] with our code. We could also broaden 
the present study using other partonic matrix elements for 
J/t// production, by considering other parametrisations for 
the shadowing. A careful comparison with results from the 
Colour-Evaporation Model at NLO [28] is planned. A bet- 
ter treatment of the interaction between the cc pair and the 
nuclear matter could be also achieved by taking into account 
coherence effects such as the c-quark shadowing [29]. Stu- 
dies of the kj broadening are also envisioned. 

In summary, we have demonstrated that the kinematics 
of the partonic processes responsible for the J/iJ/ production 
is of particular relevance to assess the importance of CNM 
effects both in pA and AA at RHIC energies. Moreover, 
we argue that a significant part of the rapidity dependence 
of Rfyji in the central collisions can be accounted by CNM 
effects only. 
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We have evaluated Cold Nuclear Matter effects on J/tfr 
production in proton-nucleus and nucleus-nucleus colli- 
sions at relativistic energies. We have considered the J/tfr 
transverse momentum effects within the specific kinematics 
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